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Abstract 

An implicit method for the ohmic dissipation is proposed. The proposed method is based on the 
Crank-Nicolson method and exhibits second-order accuracy in time and space. The proposed method 
has been implemented in the SFUMATO adaptive mesh refinement (AMR) code. The multigrid method 
on the grids of the AMR hierarchy converges the solution. The convergence is fast but depends on the 
time step, resolution, and resistivity. Test problems demonstrated that decent solutions are obtained even 
at the interface between fine and coarse grids. Moreover, the solution obtained by the proposed method 
shows good agreement with that obtained by the explicit method, which required many time steps. The 
present method reduces the number of time steps, and hence the computational costs, as compared with 
the explicit method. 
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1. Introduction 

The magnetic field plays an important role in star for- 
mation. Taking the magnetic field into account, simula- 
tions of protostellar collapse have been performed in nu- 
merous studies (reviewed by Klein et al. 2007). Most 
of these studies assumed ideal magnetohydrodynamics 
(MHD). 

Interstellar gas is partially ionized, and there are several 
processes of magnetic diffusion, e.g., the ohmic dissipa- 
tion, the Hall effect, and ambipolar diffusion. The ohmic 
dissipation is effective at high densities of n > 10 16 cm~ 3 , 
whereas the ambipolar diffusion is effective at low den- 
sities of n < 10 9 cuT 3 (e.g., Kunz & Balbus 2004). The 
timescale of the magnetic diffusion is significantly longer 
than the freefall time at n < 10 12 cm -3 (Nakano et al. 
2002), and hence magnetic diffusion does not appear to 
change the behavior of the gravitational collapse quali- 
tatively. The gravitational collapse ceases in the dense 
region of n ^ 10 11 cm~ 3 owing to the formation of an adi- 
abatic core, i.e., the first core (Larson 1969). Therefore, 
subsequently formed objects, e.g., circumstcllar disks, 
protostars, and outflows, likely suffer from magnetic dif- 
fusion. 

The recent numerical simulations for protostellar col- 
lapse begin to take into account the magnetic diffu- 
sion (e.g., Machida et al. 2006; Machida et al. 2007). 
However, the governing equation of the magnetic diffusion 
is parabolic, and therefore the time step for the magnetic 
diffusion is very small compared with the hydrodynamic 
time step when a high-resolution explicit method is em- 
ployed. High resolution is important in the simulation of 
protostellar collapse and is usually provided by means of 
adaptive mesh refinement (AMR). 



Several strategies for solving the magnetic diffusion 
have been proposed. The super time-stepping method 
is a type of explicit method, in which a large time step 
can be used (O 'Sullivan & Downes 2006; O'Sullivan & 
Downes 2007; Choi et al. 2009). However, for the diffu- 
sion dominated problem, the time step is still restricted 
to be shorter than that of the hydrodynamic time step. 
Tillcy & Balsara (2008) proposed a semi-implicit scheme 
for ambipolar diffusion using a two-fluid approximation, 
where the time step is restricted in inverse proportion to 
the drift velocity. The present author previously imple- 
mented the ohmic dissipation in a nested grid code by 
using a sub-cycle of the induction equation (Machida et 
al. 2006; Machida et al. 2007). By this method, the pro- 
tostellar collapse from a molecular cloud core to protostar 
formation was successfully simulated. Although each sub- 
cycle required a small computational cost, the number of 
sub-cycles becomes very large when solving the magnet- 
ically dissipative region, e.g., the region proximal to and 
inside of a protostar. Moreover, the resistivity was ap- 
proximated as being locally constant. 

An implicit scheme for solving the ohmic dissipation 
has been developed and implemented in the SFUMATO 
MHD- AMR code (Matsumoto 2007). In §2, the details of 
the implicit scheme are presented. In §3, the results of 
several numerical tests are presented. Finally, the paper 
is summarized in §4. 

2. Implicit scheme 

The induction equation with the ohmic dissipation is 
given by 



dB 

~dt 



Vx(oxB-j)VxB), 



(1) 
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where B, v, and 77 denote the magnetic field, velocity, and 
resistivity, respectively. Equation (1) is solved by an oper- 
ator splitting approach. The contribution of the first term 
on the right-hand side of the equation is solved explicitly 
according to Matsumoto (2007), and the contribution of 
the second term is then solved by the implicit scheme pre- 
sented herein. We hereinafter restrict our focus to the 
solution of the ohmic dissipation. 

2.1. Discretization 



and (V • f) i ■ k is calculated in the same manner. The 
differential terms in F x .i+i/2.j,k are given by 



(dxB y ) t+1 



+i/2j,fc - Ax 



(9) 



_ B x ,i+l,j+l,k + B x< i t j+i t k ~ B x ,i+l,j-l,k — £x,i,j-l,fc 



4Ay 



by 



The governing equation of the ohmic dissipation is given 



(10) 



— = -V x (77V x B) . 
dt w ; 



(2) 



The resistivity r\ at the cell surface is given by arithmetic 
average, e.g., 

Vi+l.j,k + Vi,j,k 



Vi+l/2,j,k 



(11) 



Equation (2) is written in conservation form as follows: 
r) R 

^ + V-F = 0, (3) 
where the numerical flux F = (F x ,F y ,F z ) is given by 




F x = r) 



'1 



—d x B y + d y B x 
—d x B z + d z B x 

—d y B x + d x B y 


-dyB z +d Z By 

-d z B x +d x B z 
—d z B v 







d y B z 



(4) 



(5) 



(0) 



Equation (3) is discretized as follows: 

B, jL b.^j, • AAi •: V • /*':,,/, • •:! A.A/:V •/:,,.,. 0. 

where, for convenience, the unknown variables are written 
in uppercase, and the known variables are written in low- 
ercase: B:=B n+1 ,F:=F n+1 ,b:=B n ,imdf:=F n . The 
superscript n denotes the time level, and At = t n+1 — t n . 
The subscripts i,j,k are the indexes of a cell in the x, y, 
and z directions, respectively, and are used to label cells. 
The parameter A specifies the type of temporal difference. 
The backward difference is obtained when A = 1 , and the 
central difference is obtained when A = 1/2. Therefore, 
A = 1 results in a temporal first-order accuracy, while 
A = 1/2 results in a temporal second-order accuracy. The 
case of A = 1/2 corresponds to the Crank-Nicolson scheme. 
Spatial discretization is performed with the central differ- 
ence, yielding spatial second-order accuracy. Each com- 
ponent of the numerical flux is defined at the cell surface, 
and hence the divergence of the numerical flux is calcu- 
lated as follows: 



Equation (7) is rewritten in the form of a difference 
equation as follows: 

£B%,j,k = Si,j,k> (12) 
where 

CB.^j. B, , 1. ■ \AliY-F:, !h . (13) 

S itj ,k = bi, u - (1 - A) At (V • /) . j>h . (14) 

Equation (12) indicates that the unknown B is solved by 
the linear operator C and the source term S, which is a 
function of the known b. 

2.2. Multigrid method 

Equation (12) is solved by the multigrid method. Here, 
the strategy of the multigrid method is the same as that 
of Matsumoto (2007), who solved the scalar PDE of the 
Poisson equation, while the present method solves the vec- 
or PDE of equation (2). Therefore, all of the procedures 
of Matsumoto (2007) are extended to those for vectors, 
and full- weight prolongation and averaging restriction are 
performed for each vector component. Since the smooth- 
ing procedure depends on the equation to be solved, it is 
newly developed as shown in § 2.3. 

The multigrid method £pjJj G solves R now when the ini- 
tial estimation RS UCSS and the source term S are given as 
follows: 



B r 



-C^ G (B^ SS ,S). (15) 

Since equation (12) is linear, we use the multigrid method 
itcratively, as follows: 



R=S-CB gucss 

B new = -B guess + £^ G (0 ) il) 

^ggucss ^_ ^ncw 



(V-F) 



F x .i+i/2,j,k — F xi _ 


l/2,j,fe 


Ax 




Fy,i.j + l/2,k ~ Fy,i,j 


-1/2, k 


Ay 




z,i,j,k-\-l/2 ~ -^z,i,j. 


,fc— 1/2 



Az 



(8) 



(16) 
(17) 
(18) 

This iterative utilization of the multigrid method reduces 
every component of a residual, R. 

The multigrid method given by equation (15) consists of 
(1) the full multigrid (FMG) cycle on the AMR hierarchi- 
cal grids, (2) the multilevel adaptive technique (MLAT) 
with the full approximation scheme (FAS) on these grids, 
and (3) an FMG-cycle on the base grid (Matsumoto 2007). 



No. 
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Table 1. Numbers of iterations for the multigrid method 



Schemes 


V-cycle 


Pre-smoothing 


Post-smoothing 


FMG on AMR 


4 


4 


4 


MLAT-FAS on AMR 


2 


4 


4 


FMG on base grid 


2 


2 


2 



These schemes have parameters: the numbers of itera- 
tions for V-cycle, pre-smoothing, and post-smoothing pro- 
cedures in each grid level. These parameters adopted 
through this paper are shown in table 1. These parame- 
ters affect a convergence speed of the multigrid method; 
small numbers of these iterations slow the convergence 
while the computational cost is reduced. 

As shown in section 3, several cycles of the multigrid 
method given by equations (16)-(18) reduce the residual 
by more than an order of magnitude. In the numerical 
tests, we performed 20 cycles of the multigrid method in 
order to estimate solutions converged enough. 

2.3. Smoothing 

As a smoothing operator, the red-black Gauss-Seidel 
iteration is adopted. When equation (12) is solved sepa- 
rately for Bi,j,k hi each vector component, and Bij.k is 
replaced by B"^ atcd , we obtain the following relationship: 



t = 0.0000000E+00 



^updated 



B 



i,3,k 



R 2 



■,k/{l + Oy) 

,k/(l + a z ) 



where 



R 



hj,k 




CB 



i,j,ki 



(19) 



(20) 



w f ( Vi.3-l/2M+Vi,j+l/2,k 

V Ay 2 

Vi,j,k-l/2 + Vi,j,k+l/2 

Az 2 



(21) 



Oty = XAt 



Wi,j,k-l/2 + Vi,j,k+l/2 



Az 2 

Vi-l/2,j,k + Vi+l/2,j.k 



= XAt 



Vi-l/2,j,k + Vi+l/2,j,k 



Ax 2 

Vi,j-l/2,k + r li,j+l/2,k 

Ay~ 2 



(23) 



Equation (19) gives the approximate solution of B"^ dated 
for a given initial guess of -Bj We adopt equation (19) 
as a smoothing operator. A red-black ordering is adopted 
for sweeping the grid. 




Fig. 1. Initial condition for the sinusoidal diffusion problem. 
The gray scale denotes the distribution of B z , and lines denote 
the boundaries of the AMR blocks, each of which has 32 3 cells. 



2-4- Time step 

The AMR code was equipped with two modes of time- 
marching: an adaptive and a synchronous time-step mode. 
In the former mode, a coarser grid has a longer time step 
than a finer grid, and this mode is appropriate for non- 
self-gravitational gases because the system equations are 
hyperbolic. In the latter mode, every grid-level has the 
same time step, and this mode is appropriate for self- 
gravitational gases because the Poisson equation is ellip- 
tic. For a problem including the ohmic dissipation, the 
synchronous time-step is adopted because the equation (2) 
is parabolic. 

3. Numerical tests 

3.1. Sinusoidal diffusion problem 

We consider the problem in which the sinusoidal mag- 
netic field diffuses with a constant resistivity. The initial 
magnetic field is given as follows: 



B z = sin(fc • r), 



(24) 



(22) 

and B x = B v = 0, where the wave number is set at 
k = 27r(l,2,0) T . This setting reduces the equation of the 
ohmic dissipation to the heat equation. The resistivity 
is set at r\ — 1. The computational domain is x,y,z S 
[0,1] x [0,1/2] x [0,1/4]. Periodic boundary conditions 
are imposed. The computational domain is covered by 
4x2x1 base blocks, each of which has ./V 3 cubic cells. The 
cell width is therefore given by Ax = Ay = Az = l/(4iV) 
in the base grid. The domain of x S [0, 1/2] is refined by 
blocks that are twice as fine. Figure 1 shows the initial 
distribution of B z and the block distribution for N = 32. 
The cell width in the left-hand side is Ax = 1/256, and 
that in the right-hand side is Ax = 1/128. 

We performed the convergence test by changing the 
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Fig. 2. L\ norm of error as a function of time step At for 
the sinusoidal diffusion problem. The spatial resolution is 
N = 32. The upper and lower solid lines denote the errors for 
the cases in which A = 1 and 1/2, respectively. The dashed 
lines indicate the relationships of errors in proportion to At 
and At 2 , respectively. 

time step Ai, in order to measure the temporal accuracy. 
The Li norm of the error is measured at t = 4 x 10 -3 
through comparison with the exact solution of the follow- 
ing equation: 

.B cx (r) = cxp (-ri\k\ 2 t) sin(fc ■ r). (25) 
The L\ norm is estimated as follows: 



i,j,k, 



(26) 



i,j,k 



where AVij } k denotes the volume of a cell located at r^fe, 
and V denotes the volume of the entire computational 
domain. By the stage of t = 4 x 10 -3 , the amplitude of 
B z is reduced to 0.45. 

Figure 2 shows the L\ norm as a function of the time 
step At. We examined the cases of A = 1 (backward 
difference) and A = 1/2 (Crank-Nicolson). The scheme 
with A = 1 exhibits first-order accuracy, and that with 
A = 1/2 exhibits second-order accuracy. For the scheme 
with A = 1/2, the decrease in the L\ norm with decreasing 
At is saturated at At < 2 x 10 -4 , exhibiting the constant 
L\ norm of ~ 10~ 4 . The saturation is primarily attributed 
to a discretization error. We confirmed that the value of 
the saturation decreases with decreasing cell width. 

Figure 3 shows the L\ norms estimated in the fine re- 
gion of < x < 1/2 (fine region) and the coarse region of 
1/2 < x < 1 (coarse region). For the first-order scheme 
(A = 1), the error in the coarse grid is slightly larger than 
that in the fine grid. In contrast, the second-order scheme 
(A = 1/2) exhibits an error in the coarse grid that is ap- 
proximately 6 times larger than that in the fine grid when 
A^<2xl6~ 4 . 

Figure 4 shows the distributions of the errors at At = 
10~ 4 for the schemes with A = 1 and 1/2. For the first- 
order scheme (A = 1), the error is distributed smoothly 



Fig. 3. L\ norm of error as a function of time step At for 
the sinusoidal diffusion problem. The spatial resolution is 
N = 32. The upper and lower lines denote the errors for 
the cases in which A = 1 and 1/2, respectively. The solid 
lines with diamonds and the dotted lines with filled circles 
denote the errors in x £ [0,1/2) (fine region) and x G (1/2,1] 
(coarse region), respectively. The dashed lines indicate the 
relationships of errors in proportion to At and At 2 . 



t = 4.0000000E-03 




Fig. 4. Distribution of error B z — £? e x in the x — y plane 
for the sinusoidal diffusion problem. Errors are shown for the 
N = 32 resolution and At = 10 -4 , solved by a scheme with 
A = 1 (upper) and a scheme that with A = 1/2 (lower). 

through the fine and coarse grids. For the second-order 
scheme (A = 1/2), the coarse grid shows a larger systematic 
error than the fine grid, causing the large L\ norm in the 
coarse grid. Moreover, the error is somewhat large in the 
coarse grid near the interface between the fine and coarse 
grids. 

Figure 6 shows the decrease in the maximum residual 
during the iteration of the multigrid method given by 
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Fig. 5. L± norm of error as a function of cell width 
Ax for the sinusoidal diffusion problem. Upper and lower 
solid lines indicate the errors obtained by the schemes with 
A = 1 and 1/2, respectively. The dashed lines indicate 
the relationships of errors in proportion to Af and At 2 . 

We performed a convergence test with respect to the spa- 
tial resolution by changing the number of cells inside 
a block, N 3 = 4 3 , 8 3 , 16 3 , 32 3 . The time step is set to 
At = Ax lCT 3 (4/iV) 2 . Figure 5 shows the L\ norm as a 
function of the cell width Ax for the schemes with A = 1 
and 1/2. Both schemes exhibit spatial second-order accu- 
racy. 




5 10 15 

Number of multigrid iterator! 

Fig. 6. Maximum residual |it| m ax as a function of itera- 
tion number of the multigrid method for the sinusoidal dif- 
fusion problem. Solid and dotted lines indicate data for 
A = 1/2 and 1, respectively. Diamonds and circles indicate 
data for At = 4 X 10 _ 3 and 10 -4 , respectively. Open and 
filled symbols indicate data for N = 8 and 32, respectively. 



equations (16) through (18) for various time steps At and 
spatial resolutions N. The residual is calculated according 

to equation (20), and raax(\Rx,i,j,k\A^y,i,jMA^z,i,j,k\) is 
plotted. In all cases, one iteration of the multigrid method 
reduces the residual to less than 10~ 5 . Comparison of 
the schemes with A = 1/2, and 1 reveals that the scheme 
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Fig. 7. Distribution of B z in the y = plane at t = 4 for the 
Gaussian diffusion problem with A = 1/2. Circles, diamonds, 
triangles, and asterisks denote the solutions of At = 0.5, 
1.0, 2.0, and 4.0, respectively. Green, blue, and red sym- 
bols indicate solutions on the grids of levels 0, 1, and 2, re- 
spectively. The solid curve denotes the exact solution. In 
order that all of the solutions could be plotted, the plots 
are offset from each other by 0.01 in the vertical direction. 
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Fig. 8. Same as Figure 7, but for A = 1. 



with A = 1/2 exhibits fast convergence. Moreover, there 
is a tendency whereby cases with smaller r/At/Ax 2 ex- 
hibit faster convergence. The residuals with rjAt/Ax 2 = 
262, 16.4,6.55, and 0.410 in the fine grid correspond to the 
lines with filled diamonds, open diamonds, filled circles, 
and open circles, respectively. 

3.2. Gaussian diffusion problem 

We examine the diffusion of B z in the Gaussian profile, 
the exact solution of which is given as follows: 



B z (x,y) 



1 



■ exp 



- y 



(27) 



Anr/it + to) 

where to = 1 and rj = 1. The computational domain is 
x,y,z G [—16,16], which is resolved by the base grid of 32 3 
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Fig. 9. L\ norm of error as a function of time step At 
for the Gaussian diffusion problem. Solid and dotted 
lines indicate the errors obtained by the schemes with 
A = 1/2 and 1. respectively. Diamonds, open circles, 
and filled circles denote the errors on the grids of lev- 
els 0, 1, and 2, respectively. The dashed lines indicate 
the relationships of errors in proportion to At and At 2 . 



cubic cells. The region around the z axis (x = y = 0) is 
covered by the fine grids, as shown in Figure 7. The cell 
widths are Ax = 1.0, 0.5, and 0.25 for the grids of levels 0, 
1, and 2, respectively. Periodic boundary conditions are 
imposed. 

Figure 7 shows the solutions of B z at t = 4 with var- 
ious time steps At for the scheme of A = 1/2. As the 
time step At increases, the solution deviates from the ex- 
act solution. The solution with At = 2.0 shows significant 
undulation \x\ < 4. For the solution with At = 4.0, the 
undulation mashes up the solution in the finest grid (red 
asterisks). Note that the Crank-Nicolson method is un- 
conditionally stable for the von Neumann stability analy- 
sis, while a large rjAt/Ax 2 produces such undulation due 
to violation of the maximum principle (Morton & Mayer 
2005). We also found that the undulation occurred with 
smaller At when the initial Gaussian profile had a nar- 
rower width (a smaller to). 

The scheme with A = 1 yields smooth solutions even for 
large At, as shown in Figure 8. Although monotonicity is 
maintained in the solutions, the solution in the finest grid 
deviated considerably from the exact solution when At is 
large. 

Figure 9 shows the L\ norm of the error as a function of 
time step At for A = 1/2 (solid lines) and 1 (dotted lines). 
The norm is estimated separately on each grid level. The 
errors for A — 1/2 and 1 exhibit second-order accuracy and 
first-order accuracy, respectively, on the grids of levels 1 
and 2. For the grid of level (the base grid), the de- 
pendence of the errors on At is shallower because of the 
periodic boundary conditions. Note that the scheme with 
A = 1/2 maintains second-order accuracy on the grid of 
level 2, even when considerable undulation occurs with a 
large time step. 
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Fig. 10. Distributions of the magnetic field B and the 
resistivity t] at the initial stage (left) and the stage of 
t = 1 (right), which is solved by the second-order im- 
plicit scheme with A = 1/2. Tubes and isosurfaccs 
denote the magnetic field lines and resistivity, respec- 
tively. The colors of tubes illustrate the field strength. 
The levels of the isosurfaccs are rj = 0.05, 0.425, 0.8. 
The region of x, z e [—2.4,2.4], y € [0,2.4] is shown. 



3.3. Comparison with an explicit scheme 

We compared the solutions obtained by implicit 
schemes with the solutions obtained by an explicit scheme. 
The resistivity is distributed as follows: 



r)(r) = exp [~{x 2 + y 2 + z 2 )] 
and the initial magnetic field is given by 
B z (r)^ C ^[-(x 2 + y 2 )] 



(28) 



(29) 



and B x = B y = 0. The computational domain is x,y,z S 
[—4,4], which is resolved by 64 3 cells. Periodic bound- 
ary conditions are imposed. We continue the dissipation 
process of B until t = 1 by using the implicit schemes 
with A = 1/2 and 1 and an explicit scheme. The ex- 
plicit scheme integrates time by means of the predictor- 
corrector method and a spatial central difference in order 
to achieve second-order accuracy in time and space. The 
time step is set at At = 10~ 3 for the explicit scheme, and 
At = 10 _1 for the implicit schemes. 

Figure 10 shows the initial conditions and the solution 
at t = 1 solved by the implicit scheme with A = 1/2. The 
right-hand figures show the magnetic fields bent due to 
the ohmic dissipation, indicating the reduction of B z and 
the generation of B x and B y . 

Figure 11 compares magnetic fields obtained by the ex- 
plicit scheme, the second-order implicit scheme (A = 1/2), 
and the first-order implicit scheme (A = 1). All of the 
solutions are consistent with one another. In particular, 
the solution obtained by the implicit scheme with A = 1/2 
exhibits excellent agreement with that obtained by the ex- 
plicit scheme. This is attributed to the high accuracy of 
the implicit scheme with A — 1/2, which achieves second- 
order accuracy. 

4. Summary and Discussion 

We have presented an implicit scheme for solving the 
ohmic dissipation for the SFUMATO MHD-AMR code. 
The induction equation of the ohmic dissipation is solved 




Fig. 11. Distributions of the magnetic fields at t = 1 in the y = plane. The magnetic fields are solved by the explicit scheme 
(left), the implicit scheme with A = 1/2 (middle), and the implicit scheme with A = 1 (right). Upper panels show B z , where 
the contour levels are B z = 0.05, 0.1, ■■■ ,0.95. Lower panels show B x , where the contour levels are B x = —0.09, —0.08, 0.09. 



by this implicit scheme, which is based on the Crank- 
Nicolson method, which has an option for selecting first- 
order accuracy (A = 1) and second-order accuracy (A = 
1/2) in time. For both cases, the spatially central differ- 
ence yields second order accuracy in space. 

The multigrid method is used for the convergence of 
a solution and exhibits fast convergence. Although the 
convergence speed depends on rjAt/ Ax 2 , several cycles of 
the multigrid method reduce the residual by more than 
an order magnitude. The solution is obtained over the 
AMR hierarchical grid, in which fine and coarse grids co- 
exist. Moreover, no spurious features appear at the in- 
terface between fine and coarse grids. Note that in the 
convergence process of the multigrid method, the numeri- 
cal fluxes given by equations (4) through (6) are conserved 
even at the interfaces between the fine and coarse grids, 
because of the refluxing procedure, where the fluxes of 
the coarse grid are obtained by summing those of the fine 
grids at the interface (Matsumoto 2007). This leads to 
conservation of magnetic flux. 

Since the second-order scheme of A = 1/2 is based on 
the Crank-Nicolson method, the scheme is uncondition- 
ally stable. However, the solution can contain spurious 
oscillations when rjAt/ Ax 2 is large, as shown in Figure 7. 
For example, for a one-dimensional heat equation, the 
analysis of the maximum principle leads a condition of 



rjAt/ Ax 2 < 3/2 to suppress the oscillation (Morton & 
Mayer 2005). This condition is slightly weaker than the 
CFL condition of an explicit scheme. In the astrophysi- 
cal simulations, the oscillation of the magnetic field may 
change the direction of the magnetic pressure gradient, 
which may change the phenomena of the simulations qual- 
itatively. In contrast, the first-order scheme (A = 1) retains 
monotonicity, as shown in Figure 8, but its error is larger 
than that of the second-order scheme. 

Numerical computations were carried out on the Cray 
XT4 at the Center for Computational Astrophysics, 
CfCA, of the National Astronomical Observatory of 
Japan. The present research was supported in part by 
the Hosei Society of Humanity and Environment. The 
present research was supported in part by Grants-in-Aid 
for Scientific Research (C) 20540238 and (B) 22340040 
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